Spinodal density enhancements in simulations of relativistic nuclear collisions 
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We recently introduced a fluid-dynamical model for simulating relativistic nuclear collisions in 
the presence of a first-order phase transition and made explorative studies of head-on lead-lead 
collisions. We give here a more detailed account of this novel theoretical tool and carry out more 
exhaustive studies of the phase-separation dynamics. Extracting the density enhancement caused 
by the spinodal instabilities, the associated clump size distribution, and the resulting transverse 
flow velocity, we examine the sensitivity of these quantities to the strength of the gradient term that 
promotes the phase separation, to the details of the initial density fluctuations that form the seeds 
for the subsequent amplification, and to the equation of state. 

PACS numbers: 25.75.-q, 47.75. +f, 64.75. Gh, 81.30.Dz 



I. INTRODUCTION 

The aim of high-energy nuclear physics is to elucidate 
the properties of strongly interacting matter at extreme 
temperatures and densities. A central subject in these 
investigations is the quark-gluon plasma (QGP), a new 
phase of matter in which the quarks and gluons are de- 
confined. Experiments at SPS, RHIC, and LHC have 
produced vast amounts of data indicating that such a 
new phase is indeed produced [TI4TTI| . 

Strongly interacting matter is expected to posses a rich 
phase structure. In particular, compressed baryonic mat- 
ter may exhibit a first-order phase transition that per- 
sists up to a certain critical temperature |12| . At vanish- 
ing baryon chemical potential, lattice QCD calculations 
yield a smooth transformation from confined to decon- 
fined matter at a cross-over temperature of T x w 150 
- 160 MeV [U, 0- Unfortunately, the lattice calcula- 
tions cannot be readily extended to finite chemical po- 
tential due to the so-called sign problem. While several 
methods have been employed for circumventing this basic 
problem, including Taylor expansion in fi/T, reweight- 
ing, and imaginary chemical potential, the results are not 
conclusive on the existence of a genuine first-order phase 
transition at finite chemical potential and the location of 
the associated critical point [T5l - [l8j . 

It is therefore necessary to rely on experiment for 
the determination of the phase structure of baryon-rich 
strongly interacting matter. Experimental efforts are un- 
derway to search for evidence of a first-order phase tran- 
sition and the associated critical end point |l^| - |21j ]. For 
these endeavors to be successful, it is important to iden- 
tify observable effects that may serve as signals of the 
phase structure. This is a challenging task because the 
colliding system is relatively small, non-uniform, far from 
global equilibrium, and rapidly evolving, features that 
obscure the connection between experimental observables 
and the idealized uniform equilibrium matter described 
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by the equation of state. Therefore, to understand how 
the presence of a phase transition may manifest itself 
in the experimental observables, it is necessary to carry 
out dynamical simulations of the collisions with suitable 
transport models. 

Many numerical simulations of high-energy nuclear col- 
lisions have employed ideal or viscous fluid dynamics 
which has the important advantage that the equation of 
state (EoS) appears explicitly. By contrast, in most mi- 
croscopic transport models the EoS is often unknown or 
very cumbersome to determine. The focus up to now has 
mainly been on bulk observables and their dependence 
on a softening of the EoS. For this purpose, the instabil- 
ities associated with a first-order phase transition were 
usually removed by means of a Maxwell construction, 
thereby ensuring that bulk matter remains mechanically 
stable throughout the expansion. 

However, when a first-order phase transition exists, a 
low-density confined phase (a hadronic resonance gas) 
may coexist thermodynamically with a high-density de- 
confined phase (a baryon-rich quark-gluon plasma) and, 
consequently, bulk matter prepared at intermediate den- 
sities would be unstable and seek to separate into the 
two coexisting phases. Therefore, in a nuclear collision, 
when the dynamical evolution drives the bulk density 
into the phase coexistence region, the instabilities will be 
triggered. In particular, the spinodal instabilities [22T4251 ] 
may generate a non-equilibrium evolution that in turn 
may generate observable fluctuations in the baryon den- 
sity [26l - [29| and the chiral order parameter [3(| HH . Fur- 
thermore, nucleation and bubble formation may also con- 
tribute towards the phase separation process. 

In order to ascertain the degree to which these mech- 
anisms may manifest themselves in actual nuclear col- 
lisions, we have performed numerical simulations with 
finite-density fluid dynamics, using a previously devel- 
oped two-phase equation of state and incorporating a 
gradient term in the local pressure [32| . This latter refine- 
ment emulates the finite-range effects that are essential 
for a proper description of the phase transition physics 
[22I, [2J, [25|, [33[ ■ In particular, the gradient term ensures 
that two coexisting bulk phases will develop a diffuse in- 
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terface and acquire an associated temperature-dependent 
tension. Furthermore, of key importance to the present 
study, the gradient term also causes the dispersion rela- 
tion for the collective modes in the unstable phase region 
to exhibit a maximum, as is a characteristic feature of 
spinodal decomposition [22j . Thus we employ a transport 
model that has an explicitly known two-phase equation of 
state and that treats the associated physical instabilities 
in a numerically reliable manner. 

A first application of this novel tool was recently made 
for head-on collisions of lead nuclei at various energies 
|32j . For a certain range of collision energies, several GeV 
per nucleon, the bulk of the system will reside within 
the spinodal region of the phase diagram for a sufficient 
time to allow the associated instabilities to significantly 
enhance the initial density irregularities. This charac- 
teristic phenomenon was exploited previously to obtain 
experimental evidence for the nuclear liquid-gas phase 
transition [2 21 . |3J] and we are exploring the prospects for 
it being useful for signaling the confinement phase tran- 
sition. 



II. THE EQUATION OF STATE 

In order to obtain a suitable equation of state, we em- 
ploy the method developed in Ref. (25{. Thus we work 
(at first) in the canonical framework and, for a given T, 
we obtain the free energy density fr(p) in the phase co- 
existence region by performing a suitable spline between 
two idealized systems (either a gas of pions and inter- 
acting nucleons or a bag of gluons and quarks) held at 
that temperature. In Ref. [25| the focus was restricted 
to subcritical temperatures, T < T cr it, so for each T 
the spline points were adjusted so the resulting fr(p) 
would exhibit a concave anomaly, i.e. there would be 
two densities, pi(T) and P2(T), for which the tangent of 
fr(p) would be common. This ensures phase coexistence, 
i.e. the chemical potentials match, px(pi) = Pt(P2), 
because pr(p) = dpfr(p), and so do the pressures, 
Pt(pi) = Pt(P2), because pr(p) = Pt(p)p ~ h{p)- Ref. 
[32| extended the equation of state to T > T cr it by using 
splines that are convex, as is characteristic of single-phase 
systems. After having thus constructed fr(p) for a suf- 
ficient range of T and p, we may obtain the pressure, as 
well as the energy density et(p) = fr(p) — Tdrfrip), 
by suitable interpolation and then tabulate the equation 
of state, po(e,p), on a convenient Cartesian lattice. The 
procedure is illustrated in Fig. [TJ 

Figure [5] shows a contour plot of the resulting equation 
of state in the canonical representation, pt(p)- The lower 
and upper boundaries of the phase coexistence region, 
Pi(T) and P2(T), respectively, are shown; they come to- 
gether at the critical point (/9 C rit, Tcrit)- Uniform matter is 
thermodynamically unstable inside this region, but still 
mechanically stable near the phase coexistence bound- 
aries. But inside the spinodal boundary, along which the 
speed of sound vanishes, uniform matter is both thermo- 
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FIG. 1. [Color online] The construction of the equation of 
state is illustrated for T = 0. Top pari: The free energy 
densities for either an interacting nucleon gas or quarks in 
a bag are joined by a spline, yielding a free energy density 
/t=o(p) (solid curve) having a concave anomaly, as is char- 
acteristic of a first-order phase transition. The single-phase 
partner fx(p) results when fr{p) is replaced by the common 
tangent in the coexistence density region. Bottom part: The 
corresponding plot for the pressure pr{p), which exhibits an 
undulation through the coexistence region where its Maxwell 
partner Pt (p) remains equal to the coexistence pressure. 

dynamically and mechanically unstable so, consequently, 
arbitrarily small density irregularities are amplified at a 
scale-dependent rate ^k{p,T). 

It is important to recognize that the introduction of 
the gradient term ensures that our model describes the 
instabilities not only in the unstable spinodal region, but 
also those in the surrounding metastable region in which 
finite seeds are required for amplification to occur (yield- 
ing nucleation or bubble formation [HI). Thus density 
irregularities may be amplified by the metastable region 
as well, and any clumping generated inside the spinodal 
region may be further enhanced as the system expands 
through the nucleation region. 

When the dynamical evolution is governed by ideal 
fluid dynamics (see below) the spinodal boundary is char- 
acterized by the vanishing of the isentropic sound speed, 
v s = 0, where vl = (p/h) (dp/dp) 3 , p , with s being the 
entropy density and h = p + e the enthalpy density. For 
dissipative evolutions the finite heat conductivity causes 
the instability region to widen and the boundary is then 
characterized by the vanishing of the isothermal sound 
speed, vt = 0, where v\ = (p/h) (dp/dp) T < v 2 s . These 
spinodal boundaries are also indicated on Fig. [2J 



A. Single-phase equation of state 

In order to ascertain the dynamical effects of the first- 
order phase transition, we construct a one-phase partner 
equation of state that resembles the two-phase equation 



3 




FIG. 2. [Color online] The canonical equation of state Pt(p) 
constructed by a spline procedure, as described in the text. 
The phase-coexistence boundaries (dotted lines) join at the 
critical point (solid dot). Also shown are the isothermal spin- 
odal boundary where vt = (red) and the isentropic spinodal 
boundary where v s — (blue). The familiar nuclear liquid- 
drop phase transition occurs at sub-saturation densities, while 
the confinement transition occurs at densities well above the 
saturation value p s . 



of state as much as possible. Therefore no changes are 
made outside the phase coexistence region, but for each 
subcritical temperature T < T cr ;t the two-phase free en- 
ergy density fr(p) is replaced by the value obtained by 
performing a Maxwell construction through the coexis- 
tence density region, i.e. for pi(T) < p < P2(T), namely 

f¥(p) = MPi) + (p- PiMPi) < Mp) . (i) 

We note this construction lowers the free energy density. 
This procedure is also illustrated in Fig. [TJ 

Thus, at a given subcritical temperature T, the chem- 
ical potential remains constant when the density is in- 
creased from px(T) to pi(T), because pr is the slope of 
fj^(p) which is linear. It then follows that also the pres- 
sure p = pp — f remains constant because it is linear in 
p and matches at p = Pi, fxipi) — /r(Pi) for i = 1,2. 



III. FLUID DYNAMICAL CLUMPING 

For our present investigation, we describe the evolution 
of the colliding system by ideal fluid dynamics, because 
dissipative effects are not expected to play a decisive role 
for the spinodal clumping |33| : Although the inclusion of 
viscosity generally tends to slow the growth, the dissipa- 
tive mechanisms responsible for the viscous effects also 
lead to heat conduction which has the opposite effect and 



also enlarges the unstable region (from the isentropic to 
the isothermal boundary). 

The basic equation of motion in ideal fluid dynam- 
ics expresses four-momentum conservation, d^T^ = 0, 
where the stress tensor is given by 

T^{x) = \p(x) + e{x)]u»(x)u u {x) - p(x)g^ , (2) 

where u M (x) is the four- velocity of the fluid. When taking 
account of the baryon current density, 

N»(x) = p(x)u^(x) , (3) 

the basic equation of motion is supplemented by the con- 
tinuity equation, d^N^ = 0. These equations of motion 
are solved by means of the code SHASTA [35[ in which 
the propagation in the three spatial dimensions is carried 
out consecutively. 

As mentioned above, a proper description of spinodal 
decomposition requires that finite-range effects be incor- 
porated [H, Therefore, following Refs. [25l[33|, we 
write the local pressure as 

p(r) =p ( £ (r),p(r))-a 2 e s ^V 2 ^ , (4) 

Ps Ps 

where we recall that po(e,p) is the equation of state, 
the pressure in uniform matter characterized by e and 
p. With p s = 0.153/fm 3 being the nuclear saturation 
density and e s rs m^Ps the associated energy density, 
the gradient term is normalized such that its strength is 
conveniently governed by the length parameter a. 

A. Interface tension 

As for the pressure (see Eq. ((U above), the free energy 
density is also modified by the presence of the gradient 
term. At a given global temperature T, its local value is 
always increased, 

f(r) = fr{p{r)) + We s (V^) 2 , (5) 

Ps 

where fr(p) is the free energy density in uniform matter 
at the specified density p and temperature T. Thus sharp 
changes in the density are costly and, consequently, if two 
coexisting phases of bulk matter are brought into physical 
contact a diffuse interface develops between them. This 
feature is obviously of key importance for determining the 
geometric properties of mixed-phase configurations, such 
as the preferred size of droplets or bubbles. Generally, 
the equilibrium density can be determined by minimizing 
the total free energy, subject to conservation of the total 
baryon number, 

S J lf(r) - pp(r)]dr = , (6) 

where p is the coexistence value of the chemical potential. 
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The crucial quantity is the tension associated with the 
interphase between the two coexisting bulk regions. This 
quantity is most easily determined by considering a pla- 
nar interface, in which case the above variational condi- 
tion, after insertion of the form in ([5]), yields [25[ 



d 2 
dx 2 ' 



,2 /'(•'•) 



Hr{p{x))-n 



d_ 

dp 



, (7) 



p=p(x) 



where the "depth" x is measured in the direction normal 
to the interface and pt (p) denotes the chemical potential 
in uniform matter at the specified density p and temper- 
ature T . Furthermore, in the last expression we have 
used the relation \it = d p fx to introduce the quantity 
A/(p) = fr(p) — fj?{p)i the difference between the ac- 
tual free energy density at the specified density and that 
associated with the corresponding Maxwell construction; 
it is negative for the relevant densities (those between the 
two coexistence densities). 

The solution of the above differential equation ([7]) 
yields the diffuse equilibrium density profile p(x) which 
changes smoothly from one coexistence value to the other 
as x moves through the interface region. However, due 
to the simplicity of the gradient approximation, the asso- 
ciated interface tension can readily be determined from 
the equation of state alone, without explicit knowledge 
of p(x) [25], namely 



cjt = a [2e s Af T (c)] 1/2 dc , 

J C1 (T) 



(8) 



where c = p/p s denotes the degree of compression. At 
the specified sub-critical temperature T, the integral ex- 
tends over compressions in the corresponding coexistence 
region. The interface tension ot decreases steadily as T 
is raised and finally vanishes at T cr it- It is therefore con- 
venient to characterize it by its maximum value, <jt=o ■ 

Generally the tension between two coexisting phases 
plays a key role in determining the characteristic spatial 
scale of the phase mixture that results when matter pre- 
pared in the unstable spinodal region phase separates. 
Unfortunately, the strength of the tension between nu- 
clear and quark matter is not well understood and the 
literature contains a large spread of suggested values. For 
our present studies we have adjusted the range a so that 
(To ~ lOMeV/fm 2 , which is obtained for a = 0.033 fm. 
(For comparison the surface tension of nuclear matter is 
MeV/fm 2 .) Though this value lies in the lower part 
of the range, it was preferred in a related study [28[ and 
agrees with the value of the surface tension found for the 
Polyakov-Quark-Meson model [36[ (see also the discus- 
sions in Refs. [33,111], for example). 



B. Spinodal amplification 

As mentioned above, uniform matter inside the spin- 
odal region (where v 2 < 0) is mechanically unstable and 




1.0 1.5 

Time t - t. nit [fm] 

FIG. 3. [Color online] The growth of harmonic density undu- 
lations inside the spinodal region as obtained with standard 
ideal fluid dynamics, shown by the symbols for various val- 
ues of the wave number k, together with the resulting fits to 
the expected analytical form shown by the continuous 

curves. 



density ripples of wave number k will be amplified at a 
rate jk(p> e)- Its value depends on the strength of the 
gradient term and is given by the dispersion relation, 



Ik 



,2£s P_ k A 



hpl 



(9) 



The first term reflects the familiar linear rise obtained for 
standard ideal fluid dynamics, 7^ = \o s \k, while the sec- 
ond term arises from the gradient term which introduces 
a penalty for the development of short-range undulations. 
As a consequence of these two opposing trends, the re- 
sulting growth rate jk exhibits a maximum at the favored 
length scale A op t = 27r/fc op t. 

The spinodal growth rates can be extracted by follow- 
ing the time evolution of small harmonic perturbations 
of uniform matter. Thus, imposing periodic boundary 
conditions, we consider initial systems of the form 



p(r) = p+Sp(0)s'm(kx) , e(r) 



-de(0)wa(kx) , (10) 



where (p, e) lies inside the spinodal phase region and the 
amplitudes Sp(0) and fe(0) are suitably small. Because 
the frequency is purely imaginary, u>k = H"fk, the early 
time evolution of the amplitudes will consist of growing 
and decaying exponentials having equal weights (because 
the initial state (fTU]) is prepared without any flow) pKjj . 



8p{t) « Sp(0) cosh(7 fe £) , Se(t) w Se(0) cosh( 7fc t) 



11) 



and it is then straightforward to extract the rate jk from 
the calculated amplitude growth. 

This is illustrated in Fig. [3] for the phase point (p, e) = 
(6p s ,10£ s ), which lies well inside the spinodal region, 
and using (Sp(0), Se(0j) — (0.1p s , 0.2e s ). The subsequent 
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FIG. 4. [Color online] Spinodal dispersion relation jk, the 
growth rate as a function of the wave number. The straight 
line is the exact result for gradient-free ideal fluid dynamics, 
whereas the two top curves are the corresponding numeri- 
cal results obtained with the standard grid size Aa; = 0.2 fm 
(solid) or with Ax = 0.1 fm (dashed). The four lowest curves 
were obtained with finite values of the range a as indicated. 



time evolution is obtained with ideal fluid-dynamics 
(without the gradient term for this illustration) and the 
Fourier components of the density are extracted. The re- 
sulting time-dependent amplitudes 5pk{t) are then fitted 
with the analytical form (|TT|) . As the figure brings out, 
the expected form is indeed produced, indicating that 
the numerical propagation of the unstable system is reli- 
able and that the growth rates 7^ can be extracted with 
confidence. 

Figure 2] shows dispersion relations extracted in this 
manner for various scenarios. Ideal fluid dynamics with- 
out a gradient term yields the indicated straight line, 
7fc = \vs\k, and this behavior is indeed approached by 
the numerical results when the employed spatial grid 
size Ax is reduced: Our standard value of Ax — 0.2 fm 
yields the top solid curve, while using of half that value, 
Ax = 0.1 fm, simply scales up the resulting curve by a 
factor of two, i.e. ^ki^Ax] — 27^. [Ax]. The deviation of 
the numerical 7^ from the analytical behavior \v s \k is due 
to the numerical procedure in the SHASTA code [35j . 

The other curves on Fig. 2] have been calculated with 
the gradient term included, using the indicated values of 
the range a. They all lead to curves that exhibit a max- 
imum growth rate and we see that using Ax = 0.2 fm 
ensures that the finite-a results are physically reason- 
able for range values down to a = 0.01 fm. Our stud- 
ies suggest that while the gradient-free calculations with 
Ax = 0.2 fm are reliable only up to k < 4/fm (be- 
yond which the growth rate levels off towards a con- 
stant), the use of a finite range, a > 0, extends valid- 
ity to significantly finer scales due to the suppression 
of such irregularities. Furthermore, for ranges above 
wO.Olfm the growth rates are rather insensitive to the 



employed grid size Ax, as is illustrated for our standard 
range of a — 0.033 fm (which yields a surface tension of 
(7o~10 MeV/fm 2 , as discussed above). 

In actual systems there is some degree of physical dissi- 
pation which gives rise to both viscosity and heat conduc- 
tion. To leading order, the viscosity reduces the growth 
rate by ~ ^[^r/ + C]k 2 /h, where r\ and £ are the shear and 
bulk viscosity coefficients, respectively. But the heat con- 
ductivity generally increases the growth rate (in particu- 
lar, it extends the boundary of the mechanically unstable 
region from the isentropic to the isothermal spinodal). 

IV. NUCLEAR COLLISIONS 

We have shown above that our model produces 
both meaningful interface properties, including the 
temperature-dependent strength of the associated ten- 
sion, and reasonable spinodal growth rates, including the 
emergence of an optimal phase separation scale. The 
model is therefore suitable for addressing dynamical sce- 
narios involving phase-transition instabilities. 

For a first study, we have considered central collisions 
of lead nuclei bombarded onto a stationary lead target at 
various kinetic energies, £iab- For each energy, an ensem- 
ble of several hundred separate evolutions are generated, 
each one starting from a different initial configuration 
generated by (the cascade mode of) the UrQMD model 
40 42] which treats the non-equilibrium dynamics dur- 
ing the very early stage of the collision. The switch from 
UrQMD to fluid dynamics is made at the time 

ti„it = 2i?//y£. M .-i , (12) 

where 7c. a/, is the center-of-mass frame Lorentz factor 
and R is the radius of the nucleus. At this time the 
two Lorentz-contracted nuclei have just passed through 
each other and all initial hard collisions have occurred. 
Each hadron is then represented by a gaussian of width 
Omit = 1 fm and the needed fluid-dynamical quantities, 
p(r), e(r), u(r), can then be extracted and tabulated a 
the desirable Cartesian spatial lattice [43j]. In this way, 
the effects of stopping as well as local event-by-event fluc- 
tuations are automatically included in the ensemble of 
initial fluid-dynamic states. 

A. Density clumping 

The amplification of spatial irregularities presents an 
important effect of the presence of a first-order phase 
transition [22|, [23|, [25|, [33[. In the present scenarios, 
spacial irregularities present already in the initial state, 
whereas the fluid-dynamical propagation does not gener- 
ate any spontaneous fluctuations in the course of the evo- 
lution (such fluctuations are generally produced at finite 
temperatures p4| but this refinement has not yet been in- 
corporated into the fluid-dynamical transport treatments 
of nuclear collisions). 
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FIG. 5. [Color online] Mean time evolution of the amplifi- 
cation of normalized density moments obtained for £i a b = 
3^ GeV with a = 0.033 fm (for which a w 10MeV/fm 2 ), 
using either the two-phase equation of state (solid) or its one- 
phase partner (dashed). 



FIG. 6. [Color online] Mean maximum enhancement of the 
normalized density moments for N = 7,5 (squares, circles) as 
obtained for various energies using either the two-phase equa- 
tion of state (solid) or its one-phase Maxwell partner (dashed). 



For any given event, and at any given time, parts of 
the system may lie within the unstable or metastablc 
region, and local density irregularities may then become 
amplified, whereas the rest of the matter is situated in 
a stable phase region where irregularities tend to erode. 
In order to ascertain the effect of those instabilities, we 
also carry out corresponding simulations with the one- 
phase Maxwell partner equation of state which contains 
no instabilities but is otherwise identical. 

A convenient quantitative measure of the resulting de- 
gree of "clumping" in the system is provided by the mo- 
ments of the baryon density density p(r), 

(P N ) = -Jv Jp(r) N P (r)d 3 r, (13) 

where A = Jp(r)d 3 r = (p°) is the total (net) 
baryon number. The corresponding normalized mo- 
ments, (p N ) I (p) , are dimensionless and increase with 
the order N, for a given density distribution p(r); the 
normalized moment for N = 1 is unity. 

The time evolution of the amplification of the nor- 
malized density moments is illustrated in Fig. [5] for 
N = 2, ... ,7, as obtained with both equations of state. 
The two sets of results differ strikingly: Whereas the 
one-phase equation of state hardly produces any ampli- 
fication at all, the one having a phase transition leads 
to significant enhancements, amounting to well over an 
order of magnitude for N = 7. This clearly demonstrates 
that the first-order phase transition may have a qualita- 
tive influence on the dynamical evolution of the density. 
Due to the variations in the initial conditions, the phase 
regions explored differ from one collision to the other. As 
a result, some of the evolutions experience amplifications 
that are significantly larger than the ensemble average 
(by up to a factor of five or so), whereas more evolutions 



are affected considerably less. 

Figure [5] depicts what would happen if the expansion 
were to continue while maintaining local equilibrium so 
the generated density enhancements eventually subside. 
However, as the hadronic gas grows more dilute, local 
equilibrium cannot be maintained. Consequently, if the 
decoupling occurs sufficiently soon after the clumps are 
formed, the associated phase-space correlations may sur- 
vive. This issue requires more detailed studies which are 
underway. 

The degree of density clumping generated during a col- 
lision depends on how long time the bulk of the matter 
is exposed to the spinodal instabilities. The optimal sit- 
uation occurs for collision energies that produce maxi- 
mum bulk compressions lying well inside the unstable 
phase region because the instabilities may then act for 
the longest time [25|, |H, [H| . At lower energies an ever 
smaller part of the system reaches instability and the re- 
sulting enhancements are smaller. Conversely, at higher 
energies the maximum compression occurs beyond the 
spinodal phase region and the system is exposed to the 
instabilities only during a relatively brief period during 
the subsequent expansion. For still higher energies the 
spinodal region is being missed entirely. 

Figure [6] shows the (ensemble average) maximum en- 
hancement achieved as a function of the beam energy 
for the two equations of state. The existence of an op- 
timal collision energy is clearly brought out. While the 
presently employed equation of state suggests that this 
optimal range is E\ a \, s» 2 — 4 A GeV, it should be rec- 
ognized that others may lead to different results. Fur- 
thermore, the magnitude of the effect depends on the de- 
gree of fluctuation in the initial conditions which in turn 
are governed by the pre-equilibrium dynamics. On the 
other hand, our studies suggest that the optimal energy 
is rather insensitive to the range parameter a. 
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FIG. 7. The size distribution of the density clumps produced 
in central lead- lead collisions at 3 A GeV using p m m = 7p s to 
define the clump boundary. The solid histogram shows the 
distribution obtained for the two-phase equation of state and 
the solid line represents an exponential fit. The distribution 
obtained with the one-phase equation of state is shown by 
the dotted histogram and the difference between the two size 
distributions is also depicted. 
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FIG. 8. (Color online) The time evolution of the average 
transverse flow velocity U±_ for central lead-lead collisions at 
a fixed-target beam kinetic energy of E\ a h =3A GeV for the 
two-phase equation of state (solid) and its one-phase Maxwell 
partner (dashed). The arrows indicate the average completion 
of the fluid-dynamical stage, defined as the time when all 
fluid cells have energy densities below five times the nuclear 
ground-state energy density. 



B. Size distribution 

To gain a more detailed understanding of the clump- 
ing phenomenon, we have studied the distribution of the 
clump sizes. Although the "clumps" tend to remain fairly 
diffuse, we may define their extension by means of a spec- 
ified density cutoff, p m i n , and then extract the total net 
baryon number contained within the resulting volume. 
Figure [7] shows the size distribution obtained for a den- 
sity cutoff of yOmin = 7p s , for central lead- lead collisions 
at 3^ GeV. 

The initial size distribution is approximately exponen- 
tial and that feature is well preserved during the evo- 
lution with the one-phase equation of state which pro- 
duces negligible amplification. The spinodal instabilities 
in the two-phaseequation of state leads to a preferential 
amplification of length scales near the optimum size, as 
is brought out by the difference between the two-phase 
size-distribution and the one obtained with the Maxwell 
partner; this difference peaks at clumps containing 5-8 
baryons. Nevertheless, for a wide intermediate range, 
from about 3 to about 16, the resulting two-phase size 
distribution retains an approximately exponential ap- 
pearance, but with a significantly gentler slope. 



C. Transverse Flow 

Radial flow has long been thought to be a affected 
the presence of a first-order phase transition because of 
the accompanying softening of the equation of state (45[ . 
It may therefore be expected that the transverse flow 



exhibits a corresponding weakening. However, recent in- 
vestigations using state-of-the-art dynamical models sug- 
gest that a simple connection between the softening of 
the equation of state and the radial flow (and its angular 
moments) cannot so easily be drawn (4(| [4tJ . Although a 
softer equation of state does reduce the build up of flow 
, it also leads to a longer expansion time which tends to 
have an opposite effect. 

We investigate here the effect of the phase transition on 
the transverse flow to determine whether it may serve as 
a useful signal. For this purpose, we define the transverse 
flow velocity U± as 



U±(t) = 



1 



M{t)Q{t) 



[v x (r, t)x + v y (r, t) y] p(r, t) d 3 r, 

(14) 



where p(r) is the local (net) baryon density at the po- 
sition r = (x,y, z), while v x (r) and v y (r) are the trans- 
verse components of the local fluid velocity v(r). M and 
Q are defined as follows, 



M(t) 
M(t)Q(tf 



p(r,t)d 3 r, (15) 
[x 2 + y 2 ]p{r,t)d 3 r . (16) 



Figure [8] shows our results for the average transverse 
flow of baryons as a function of time for central lead- 
lead collisions at the fixed-target beam energy of -Ei a b = 
3 A GeV for which the density moments exhibited their 
largest enhancement. Indeed transverse flow builds up 
more slowly for the two-phase equation of state than for 
its one-phase Maxwell partner, though the difference be- 
tween the two cases is very small, probably too small to 
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a = 0.05 fm 




Baryon number in cluster 

FIG. 9. (Color online) The size distribution of the density 
clumps produced in central lead-lead collisions at 3 A GeV as 
obtained for various values of the range parameter a with the 
two-phase equation of state. The dashed line shows an expo- 
nential fit to the standard result (obtained for a = 0.033 fm) 



be of experimental utility. Furthermore, the larger stiff- 
ness of the one-phase equation of state causes the expan- 
sion to be faster, so the end of the the fluid-dynamical 
stage (the time when the energy density has become suf- 
ficiently low, e = 5e s in the present case) will occur 
sooner and this shorter expansion time will largely off- 
set the (small) gain in flow velocity. As a result, the av- 
erage transverse flow velocities become almost identical 
at the end of the fluid-dynamical stage. Consequently, 
this quantity does not seem promising as a signal for the 
phase transition. 



V. PARAMETER STUDIES 

The fluid dynamical model employed has several pa- 
rameters that are not well constrained by our current 
knowledge of QCD. These include the strength of the 
interface tension, the spatial scale of the initial density 
irregularities, and the equation of state of baryon-rich 
uniform matter, In the following, we investigate the de- 
pendence of our results on these quantities by varying 
the respective controlling parameters one at a time. 



A. Finite range 

Figure [9] shows the dependence of the clump size dis- 
tribution on the range parameter a governing the inter- 
face tension between baryon-rich confined and deconfincd 
matter, for the optimal collision energy. We have shown 
above that density irregularities in the unstable phase re- 
gion will grow faster for smaller ranges. This feature is 
also reflected in the size distribution: a smaller a value 
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FIG. 10. (Color online) The time evolution of the seventh 
normalized density moment for central lead-lead collisions at 
E\ a b ~ 3A GeV for three different values of the range param- 
eter a, using the standard value atmt = 1 fm. 

enhances the yield of larger clumps (i.e. density enhance- 
ments with larger baryon contents) and, conversely, a 
larger a value decreases the large- A yield. 

Figure [10] shows how the variation in the range param- 
eter a affects the time dependence of the seventh normal- 
ized density moment. As expected, a smaller value of the 
range a leads to a stronger enhancement of irregularities, 
as already brought out in figure [9] 

B. Initial fluctuations 

Figure [UJ shows the dependence of the size distribution 
on the parameter CT; n it characterizing the spatial scale of 
the initial density irregularities. Although a reduction of 
the scale from the standard value of <r init = 1 fm shifts 
the clump sizes to larger values, the difference between 
the cases of <7i n it = 0.7 fm and <r- m it = 0.5 fm is very 
small, indicating that the fluctuation growth saturates 
when the scale of the initial irregularities is reduced. 

Figure Q2] shows how the reduction of <7j n i t affects the 
evolution of the seventh normalized density moment. 
The behavior is somewhat complicated. While the ini- 
tial value of the moment increases when the fluctuation 
scale (Ti n it is reduced, its relative growth tends to be- 
come smaller. This feature may be ascribed to the fact 
that a system with larger initial density fluctuations also 
has larger density gradients, which in turn implies that 
also a smaller part of the system lies inside the unsta- 
ble phase region. Furthermore, the larger gradients also 
cause the expansion through the unstable coexistence re- 
gion to proceed more rapidly. 

For our present studies, the initial state of the fluid- 
dynamical evolution has been obtained by use of the 
UrQMD model in its cascade mode in which there are 
neither inter-nuclear potentials nor nuclear mean fields. 
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FIG. 11. (Color online) The size distribution of the clumps 
produced in central lead-lead collisions at 3 A GeV obtained 
for different values of the parameter crinit controlling the spa- 
tial scale of the initial irregularities, using a — 0.033 fm and 
hte two-phase equation of state. The dashed line shows an ex- 
ponential fit to the standard result (obtained for aiuit = 1 fm). 



However, for nuclear collisions at the relatively low en- 
ergies of relevance here, below £J lab s» 5 — 10 A GeV, a 
realistic description of the nuclear interactions is impor- 
tant for describing observables like the elliptic or directed 
flow (48|. An UrQMD treatment of the initial state that 
includes mean-field type nuclear interactions [49j will de- 
crease the initial compression of the fireball. 

Figure [13] compares the maximum enhancements of 
the fifth and seventh normalized density moments ob- 
tained with the cascade-mode initial conditions used in 
the present study with those obtained when nuclear inter- 
actions are included. It is obvious that an enhancement 
is obtained only when the initial density has been gen- 
erated in the cascade mode. This can be understood as 
a result of the considerably decreased initial compression 
due to the nuclear repulsive force which prevents the bulk 
of the system achieving the compressions characteristic 
of the unstable phase region. 

The results bring out the fact that the degree of den- 
sity enhancement caused by the phase transition depend 
strongly on the initial compression relative to the loca- 
tion of the unstable region in the phase diagram. 
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Elapsed time t - t init [fm] 

FIG. 12. (Color online) The time evolution of the seventh 
normalized moment for central lead-lead collisions at 3A GeV 
for three different values of the parameter erinit controlling the 
scale of the initial density fluctuations, using a = 0.033 fm. 



It is then readily shown that the associated temperature 
and chemical potential are given by 



T b c(s,p) = -T (be,cp) , p bc (e,p) 



-p {be,cp) , (17) 



where To(e,p) — l/d e oo(s,p) is the temperature of 
the standard equation of state, while the corresponding 
chemical potential is po(e,p) = —Ta{e 1 p)/d p ao{e 1 p). It 
follows that the pressure p = To — e + pp is given by 



Pbc(e,p) = -p (bs,cp) . 



(18) 



The enthalpy and the free energy scale similarly, i.e. 
h bc (e,p) = h (be,cp)/b and fbc{s,p) = fo(be,cp)/b. It 
is also easy to see that the compressibilities remain in- 
variant, c b v c (e : p) = c"(6e,cp) and c bc (e,p) = Cp(be,cp). 
The same holds for the sound speeds, 

v b s c (e,p) = v° s (be,cp) , v b T c (e,p) = v° T (be,cp) , (19) 

thus ensuring that the isentropic and isothermal spinodal 
boundaries stretch in proportion to the employed scale 
factors. In particular, the scaled critical point is given 

by (4f,PtT) = (^ it /b, P f t /c). 



C. Equation of State 

To further illuminate this latter feature, we now modify 
the equation of state and investigate the dependence of 
the results on the location of the unstable region. 

For this purpose, we consider a family of equations of 
state obtained by rescaling the basic mechanical densities 
e and p. Thus, a particular scaled equation of state has 
the entropy density <Tb c (e,p) = ero(be,cp), where oo(e,p) 
is the entropy density of our standard equation of state. 



The scaling properties exhibited above provides a con- 
venient means for exploring different equations of state. 
As an illustration, Fig. [TJ] compares the energy depen- 
dence of the density enhancements obtained with our 
standard equation of state to those obtained with an 
equation of state that has been stretched upwards by 
30% in both directions, i.e. b = 1/1.3 and c = 1/1.3. 
We see that the qualitative features remain very similar. 
In particular, there is still an optimal collision energy. 
The larger enhancement obtained for the scaled equation 
of state is primarily due to the corresponding expansion 
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FIG. 13. (Color online) The collision energy dependence of 
the average maximum enhancement of the normalized density 
moments for N = 7, 5 as obtained for initial states generated 
either by the cascade mode of the UrQMD model (our stan- 
dard) or with the inclusion of mean fields. 
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FIG. 14. (Color online) The average maximum enhancement 
of the normalized density moments for N = 7, 5 as a function 
of the beam kinetic energy Ei a b, obtained with either our 
standard equation of state (solid curves) or the scaled one 
corresponding to b = 1/1.3 and c = 1/1.3 (dashed curves). 



of the unstable phase region which allows the amplifica- 
tion to proceed for a longer time. However, the optimal 
collision energy is only slightly increased. 



VI. SUMMARY 

As reported recently [Hj], we have augmented an ex- 
isting finite-density ideal fluid dynamics code with a gra- 
dient term and thereby obtained a transport model that 
is suitable for simulating nuclear collisions in the pres- 
ence of a first-order phase transition. It describes both 
the temperature-dependent tension between coexisting 
phases and the amplification of the spinodal modes. Ap- 
plying this novel model to lead-lead collisions, using an 
equation of state with a first-order phase transition, we 
found that the associated instabilities may cause signifi- 
cant amplification of initial density irregularities, relative 
to what would be obtained without the phase transition. 
Because such clumping may give rise to observable ef- 
fects that could signal the presence of the phase transi- 
tion, perhaps through enhanced composite-particle pro- 
duction, it is interesting to explore this phenomenon in 
more detail. 

We have here given a more detailed account of this 
novel simulation tool and carried out more exhaustive 
studies of the phase-separation dynamics. In particular, 
we extracted the density enhancement, the clump size 
distribution, and the transverse flow velocity and exam- 
ined the sensitivity of these quantities to the strength of 
the gradient term that promotes the phase separation, 
to the details of the initial density fluctuations that form 
the seeds for the subsequent amplification, and to the 
equation of state, all quantities that are yet not very well 
understood. We found that our results are robust against 
moderate changes of these model parameters. 

While the presence of a phase transition greatly in- 
creases the density enhancements, it has only little effect 
on the resulting transverse flow. Perhaps most impor- 
tantly, studies using different equations of state support 
the general existence of an optimal collision energy range 
within which the phase-transition instabilities have the 
largest effects on the dynamical evolution. Our results 
suggest that this energy corresponds to several GeV per 
nucleon of kinetic energy for a fixed-target configuration, 
a range that may be too low to access effectively at RHIC 
but which should match well with both FAIR and, espe- 
cially, NICA. 
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